Universal Dynamics of Phase-Field Models for Dendritic Growth 
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We compare time-dependent solutions of different phase-field models for dendritic solidification 
in two dimensions, including a thermodynamically consistent model and several ad hoc models. 
The results are identical when the phase-field equations are operating in their appropriate sharp 
interface limit. The long time steady state results are all in agreement with solvability theory. No 
computational advantage accrues from using a thermodynamically consistent phase-field model. 
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Dendrites are the most commonly observed solidifica- 
tion microstructures in metals. The free growth of a 
single dendrite is a prototype for problems of pattern 
selection in materials science jl|-|5| and has been exten- 
sively studied experimentally and theoretically. It is still 
not possible to compare theory with experiment because 
of the difficulties in computing three dimensional mi- 
crostructures in the temperature and material's parame- 
ter range of the experiments. 

Recently, a significant step forward was taken by 
Karma and Rappel |||7j], who not only showed how to 
compute accurately two dimensional dendritic growth 
but also were able to compare their results with theo- 
retical predictions. Their calculations used the so-called 
phase-field formulation of solidification, in which a math- 
ematically sharp solid-liquid interface is smeared out or 
regularized and treated as a boundary layer, with its own 
equation of motion. The resulting formulation, described 
in more detail below, no longer requires front tracking 
and the imposition of boundary conditions, but must be 
related to the sharp interface model by an asymptotic 
analysis. In fact, there are many ways to prescribe a 
smoothing and dynamics of a sharp interface consistent 
with the original sharp interface model, so that there is 
no unique phase-field model, but rather a family of re- 
lated models. An important part of Karma and Rappel's 
work was an improved asymptotic analysis which allows 
coarser spatial grids to be used in the numerical compu- 
tations than was previously possible. 

Although the phase-field method has gained accep- 
tance as a useful way to study solidification problems, 
a debate still exists over the interpretation and validity 
of the phase-field models themselves. Each model in- 
cludes a double-well potential field which enforces the 
above properties of the phase-field. Some models can 
be shown rigorously to satisfy entropy inequality H|). 
These are sometimes called "thermodynamically consis- 
tent" models. On the other hand, it has been argued 
that the precise form of the phase-field equations should 
be irrelevant so long as the computations are performed 



at the asymptotic limit where the phase-field model con- 
verges to the sharp interface limit jL0|. 

The purpose of this Letter is to compare the dynamics 
of the different phase-field models proposed. To this end, 
we have performed accurate and extensive computations 
using a specially developed adaptive mesh refinement al- 
gorithm [ pdp^ ]. We find that when properly used, all 
phase-field models give precisely equivalent results; not 
only does each phase-field model converge to the steady 
state predicted by theory, but also the transient dynam- 
ics approach the steady state uniquely. Indeed, once one 
has established that there is genuine universal dynamic 
behavior, the only remaining consideration is the compu- 
tational efficiency. Our results clearly indicate that the 
CPU times required for the different models are identical. 
In particular, we find no advantage for the thermodynam- 
ically consistent model. 

A secondary purpose of this Letter is to detect the 
limit of the validity of phase-field models in describing 
the sharp interface problem. In the context of the asymp- 
totic analysis of Karma and Rappel J6|,(7| , the ratio of the 
interface width to the diffusion length (referred to as the 
interface Peclet number, IPe) must be small in order for 
the different phase field-models to collapse to identical 
free-boundary problems. We show how finite-IPe dis- 
crepancies encountered between different models can be 
eliminated by adjusting the phase-field parameters. We 
emphasize that IPe is a free parameter and can be var- 
ied for numerical convenience by changing the interface 
width. 

The solidification of a pure substance is described by 
a free-boundary problem for the temperature in the solid 
and liquid phases, and the position of the interface be- 
tween them: 

d t u = DV 2 u (1) 
V n = D(d n u\+ - d n u\-) (2) 
Ui = -d(n)K - f3(ri)V n (3) 

The temperature T has been rescaled as a dimensionless 
thermal field u = (T — T m )/(L/C p ), where T m , L, and 
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Cp represent the melting temperature, the latent heat 
of fusion, and the specific heat at constant pressure, re- 
spectively. The thermal diffusivity D in Eq. (1) is as- 
sumed to be equal in both phases. Eq. (2) describes 
energy conservation at the solid-liquid interface, where 
V n is the local outward normal interface velocity and d n 
refers to the outward normal derivative at the interface 
for the solid (+) and liquid (-) phases. Finally, Eq. (3) is 
known as the Gibbs-Thomson condition, describing the 
deviation of the interface temperature Ui from equilib- 
rium, due to the local curvature k, and interface kinetics. 
d(ft) = ^{n)T m C P l ' L? is the anisotropic capillary length, 
proportional to the surface tension 7(7?), and (3 (ft) is the 
anisotropic kinetic coefficient. 

Eqs. (l)-(3) have been studied extensively to determine 
the steady state features of dendritic growth . These 
equations admit a family of discrete solutions. Only the 
fastest growing of this set of solutions is stable, and this 
is the dynamically selected "operating state" for the den- 
drite, corresponding to a unique tip shape and tip veloc- 
ity. This theoretical treatment is usually called solvabil- 
ity theory. Recent calculations of dendritic growth using 
phase-field models have been found to be in good agree- 
ment with the predictions of solvability theory P,[7pl|]. 

The phase-field model finesses the computational diffi- 
culties associated with front-tracking on a discrete lattice 
by introducing an auxiliary continuous order parameter, 
or phase-field, </>(r, t) that couples to the evolution of the 
thermal field. The dynamics of 4>(r, t) are designed to fol- 
low the evolving solidification front |p^|-|l9|. The phase- 
field interpolates between the solid and liquid phases, at- 
taining a different constant value in each phase (typically 
±1), with a rapid transition region in the vicinity of the 
solidification front. The liquid-solid interface is defined 
by the level set of </>(r, t) — 0. 

We consider phase-field equations of the form 



(4) 



(5) 



as in Refs. ||[7|. The order parameter is defined by 0, 
with <fi — +1 in the solid, and <fi = — 1 in the liquid. The 
interface is defined by <j) = 0. 

The function F(<j>, Xu) = f((f>) + Xug((f>) is a phe- 
nomenological free energy where f(cf>) has the form of 
a double-well potential, A controls the coupling between 
u and <p, and the relative height of the free energy min- 
ima is determined by u and g(4>). The function h(cj>) ac- 
counts for the liberation of latent heat. Anisotropy has 
been introduced in Eq.(5) by defining W(ft) = W A(ft) 
and r(n) = t A (ft). t q is a time characterizing atomic 



movement in the solid-liquid interface, W is a length 
characterizing the width of the interface, and 



A(ft) = (1 - 3e) 
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with A(ft) e [0, 1]. <j). x and <p <y represent partial deriva- 
tives with respect to x and y, and the vector ft = 
(4>, x x + 4> ,y%l) I {4> 2 X + ^y) 1 ^ 2 is the normal to the contours 
of <f). The constant e parameterizes the deviation of W(ft) 
from W and is a measure of the anisotropy strength. 

We use the asymptotic relationships of Karma and 
Rappel HQ to map the phase-field model into the free- 
boundary problem, where Eqs. (4) and (5) reduce to 
Eqs. (l)-(3). In terms of A(ft), (3 (ft) = f3 A(ft) and 
d(ft) — d Q [A(ft) + dgA(ftj\ , where 6 is the angle be- 
tween ft and the a;-axis; noting that tan(#) = 4> t y/4> tX , 
these expressions become (3 (ft) = f3 (l + e cos AO) and 
d(ft) = d Q (l — 15ecos46>) in the free-boundary prob- 
lem. The parameters of the phase-field model are re- 
lated to the free-boundary parameters by A = W a,i/d 
and t = W^aia2/(d D D) + W%(3 /d - The positive con- 
stants a\ and a-i depend on the exact form of the phase- 
field equations. In choosing to simulate particular mate- 
rial characteristics, we fix the experimentally measurable 
quantities d Q , f3 a , and D, leaving W as a free parameter 
which determines A and r Q . 

We compute four-fold symmetric dendrites in a 
quarter-infinite space using a new finite-element adap- 
tive grid method reported in Refs. | pT|,p!2[ . Solidification 
is initiated by a small quarter disk of radius R centered 
at the origin. The order parameter is initially set to its 
equilibrium value (/> (x) = — tanh((|x| — ii )/v2) along 
the interface. The initial temperature is u — in the solid 
and decays exponentially from u = at the interface to 
u = — Aasaf— >oo, where the far-field undercooling is 
A = (T m — Too) / (L / C p ) and T^ is the temperature far 
ahead of the solidification front in the liquid. 

The different phase-field models we study are summa- 
rized in Table B. To satisfy the asymptotics, /(</>) is cho- 
sen to be an even function, and g(4>) and h(<p) are odd. 
For all of the models, f((f>) = </> 4 /4 — </> 2 /2. For com- 
putational purposes, the g(cj>) are chosen such that the 
two minima of F(<p,Xu) are fixed at <fr = ±1. Model 1 
is a form used by AlmgrenpCj, Model 2 is a form used 
by Karma and Rappel llf^^and Model 5 is the ther- 
modynamically consistent form used by Wang et al ||. 
Models 3 and 4 are forms created by us that meet the 

TABLE I. Summary of phase-field models studied. 
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TABLE II. Summary of simulation parameters. 
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2 


1.5 


0.0469 


0.031 



above requirements. It should be noted that Model 1 re- 
quires that A < 1 /A otherwise the cf> — —1 state becomes 
linearly unstable. 

In our simulations, the computational domain is an 
L x L square box. Computations were performed at 
A = 0.65, 0.55, and 0.45. A summary of the parame- 
ters used for each simulation is given in Table ||, where 
V = Vd /D is the dimensionless tip velocity predicted 
by solvability theory, Ax is the minimum grid spacing 
of our mesh |ll],[l^], and At is the simulation time step. 
The phase-field parameters were chosen for each model so 
that they all simulated the same free-boundary problem. 
For all simulations e = 0.05, (3 = 0, and W = 1. Fig- 
ure 1 shows the dimensionless tip velocity of the dendrite 
versus time for the simulations performed at A = 0.55 
and 0.45. These results show that all of the phase-field 
models studied produce identical results for the entire 
temporal evolution of the dendrite and also converge to 
steady state solutions that are within a few percent of 
those predicted by solvability theory. In addition, the 
CPU times required for each of the models were identi- 
cal. 
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FIG. 1. Time evolution of the dimensionless tip velocity for 
five different phase-field models at A = 0.55 and 0.45. Each 
curve consists of five solutions superimposed on one another. 

At A = 0.65 (with d Q = 0.5), however, there are signif- 
icant quantitative differences between the various phase- 
field models, as shown in Fig. 2. This discrepancy is 
attributed to finite-IPe corrections at higher order in the 
asymptotic expansion. The deviations are a signal that 



the solutions have not converged as a function of the ex- 
pansion parameter and that the phase-field equations are 
not operating within the sharp interface limit. It should 
be possible to make these higher order terms negligible 
if one makes IPe smaller. Figure 2 also shows that each 
model has different convergence properties. However, in 
other simulations we have found that no single model 
consistently converges more rapidly than the others; in 
general, the convergence appears to depend on the initial 
conditions. 
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FIG. 2. Time evolution of the dimensionless tip velocity 
for five phase-field models at A = 0.65 with d = 0.5 and 
IPe= 0.094. When IPe is too large the different models do 
not exhibit universal behavior. 

In their asymptotic analysis, Karma and Rappel ex- 
pand the solutions to the phase-field equations treat- 
ing the interface Peclet number as a small parame- 
ter ^U^. This expansion parameter can be written as 
IPe = WV/D — W V/d , where the dimensionless ve- 
locity V depends only on A and e. We decreased IPe for 
the A = 0.65 case by simulating with a larger value of d a . 
Figure 3 shows the tip velocity versus time for the dif- 
ferent phase-field models at A = 0.65 with the capillary 
length d Q = 1.5. The universal behavior of the different 
models is recovered after this adjustment is made. 

We have demonstrated that one can obtain identical 
results from different phase-field models by choosing the 
expansion parameter IPe to be sufficiently small. Un- 
fortunately, in practice, the interface width is the only 
parameter that can be used to control the size of IPe, 
since V is fixed for a given A and e, and d is set by 
the particular material to be simulated. Thus, there is 
only the one free parameter, W , that can be adjusted 
make IPe smaller. This restriction can hinder computa- 
tional efficiency, as the number of grid points necessary 
to resolve the interface (and thus the simulation time) 
scales as 1/W% on an adaptive grid, and as 1/Wg on 
a fixed grid. In addition, with zero interface kinetics, 
t ~ W% which places a restriction on the computational 
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FIG. 3. Time evolution of the dimensionless tip velocity 
for five phase-field models at A = 0.65 with d = 1.5 and 
IPe= 0.031. Decreasing IPe recovers universal behavior. 

time step if an explicit scheme is used. For the simula- 
tion at A = 0.65 with d a = 0.5, an IPe = 0.031 can be 
obtained by reducing W Q by a factor of 3, but this would 
require an impractical amount of computing time. We 
note that the asymptotics of Karma and Rappel become 
most accurate at lower undercoolings fl2|| , which is also 
an experimentally relevant regime. At low A, V ~ A 4 , 
allowing the use of larger values for W /d . Simulating at 
low A requires larger system sizes. This, however, adds 
very little computational complexity for adaptive mesh 
based codes. 

We can extend the range of validity for these phase- 
field models by carrying out the asymptotic analysis fur- 
ther so that finite-IPe corrections are pushed to higher 
orders. This will lessen the restrictions on the interface 
width, thus rendering the phase-field approach compu- 
tationally more efficient. Detailed results will be pre- 
sented in a forthcoming paper pl[ . There appears to 
be a general trend that as one goes to higher orders in 
the asymptotic expansion more constraints are required 
on the functions f(4>),g(<p), and h((p) in order to get rid 
of correction terms inconsistent with the free-boundary 
formulation. These constraints can cause the phase-field 
to have solutions that are not monotonic in the inter- 
facial region, thus requiring higher grid resolution and 
computation time ^0|. We are currently pursuing the 
development of a phase-field model from a renormaliza- 
tion group approach with the goal of creating a more 
systematic convergence to the free-boundary problem. 
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